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Abstract 


High order horizontal diffusion of the form KV 2m is widely used in spectral models 
as a means of preventing energy accumulation at the shortest resolved scales. In the 
spectral context, an implicit formulation of such diffusion is trivial to implement. The 
present note describes an efficient method of implementing implicit high order diffusion 
in global finite difference models. The method expresses the high order diffusion equa- 
tion as a sequence of equations involving V 2 . The solution is obtained by combining 
fast Fourier transforms in longitude with a finite difference solver for the second order 
ordinary differential equation in latitude. 

The implicit diffusion routine is suitable for use in any finite difference global model 
that uses a regular latitude/longitude grid. The absence of a restriction on the timestep 
makes it particularly suitable for use in semi-Lagrangian models. The scale selectivity 
of the high order diffusion gives it an advantage over the uncentering method that has 
been used to control computational noise in two-time-level semi-Lagrangian models. 
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1 Introduction 


Diffusion and filtering have been common features of numerical models since the beginning 
of numerical weather prediction (Shuman, 1957). The reason for this is two-fold. First, 
diffusion and filtering can serve as crude parameterizations of subgrid scale processes that 
are not represented by the model. In the real atmospheric flow, there is a cascade of energy 
to smaller and smaller scales due to nonlinear interactions, until the motion is dissipated 
by friction. But in a model, even with the largest computer available, only scales far larger 
than the inertial subrange of turbulence can be resolved. In other words, what needs to 
be parameterized is much more than molecular diffusion. Secondly, models are usually 
subject to undesirable high wavenumber noise due, for example, to numerical discretization 
procedures. Diffusion and filtering can serve to reduce or eliminate this kind of noise. 

Explicit formulations of diffusion equations are only conditionally stable (Mesinger and 
Arakawa, 1976). As a consequence, the coefficient of diffusion cannot be increased beyond a 
critical value that depends on gridlength and timestep. This imposes a stringent restriction 
in application to atmospheric models, especially those such as semi-Lagrangian models that 
allow the use of large timesteps: the allowed diffusion coefficient is too small to eliminate 
the high wavenumber noise. Such a restriction is particularly marked in global models using 
a uniform latitude/longitude grid where the convergence of the meridians leads to a grid 
size much smaller at high than at low latitudes. A diffusion coefficient that is barely large 
enough to ensure smooth solutions at low latitudes can far exceed the critical value for 
computational stability at high latitudes. This restriction becomes even more marked if a 
diffusion equation of higher order is used. 

An alternative to diffusion are spatial filters; see Raymond and Garder (1991) for a review. 
The most commonly used in meteorological models is the Shapiro filter (Shapiro, 1971), 
which has the advantages of being very simple in formulation and of being easily extendible 
to higher orders. The filter is designed in a manner analogous to that used in solving 
the explicit formulation of a diffusion equation of equivalent order with the computational 
stability criterion satisfied. When large timesteps are used, the implied diffusion coefficient 
is so small that high wavenumber noise cannot be effectively reduced. In addition, the 
preferred spherical formulation of the filter, which involves simply averaging among the 
surrounding points, does not preserve area- weighted mean values. This implies violation of 
global conservation of the quantity being filtered. All attempts at using the filter in the 
context of spherical geometry have been shown to lead to degraded performance at high 
latitudes (Shapiro, 1979). 

In some semi-Lagrangian models, an uncentered discretization has been employed to damp 
computational noise (Tanguay, ei a/., 1992; Bates, et a/., 1993; Gravel, et a/., 1993). This, 
however, can reduce the accuracy of those terms discretized semi-implicitly (the second 
order accuracy is decreased to first order). Also, being insufficiently scale selective, the 
uncentering can have a negative impact on the large scale features in both amplitude and 
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phase, a s has been recognized in our experience and also by Semazzi and Dekker (1994). 

Implicit high order horizontal diffusion is trivial to implement in spectral models since the 
model variables are expressed in spectral harmonics which are eigenfunctions of the spherical 
Laplace operator. In this context, implicit diffusion is widely used (e.g., Williamson and Ol- 
son, 1994; Ritchie, et a/., 1994). The need for an implicit formulation of horizontal diffusion 
in gridpoint models has also been recognized, but the solution is much less straightforward. 
Jakimow et al. (1992) presented an implicit formulation for second order horizontal diffu- 
sion, and showed that its application to a semi-Lagrangian hemispheric model gave positive 
results. It was found that the optimum value for the diffusion coefficient exceeds the value 
allowed by the stability criterion of an explicit formulation. McDonald (1994) has devel- 
oped an implicit solver for an approximated form of V 2m diffusion for limited area finite 
difference models, and has found the diffusion to be effective in reducing noise in the semi- 
Lagrangian context. Our experience with the NASA/GLA Semi-Lagrangian GCM (Bates 
et al. } 1993; Moorthi et a/., 1994) indicates that at moderate spatial resolution the intrinsic 
smoothing associated with the interpolations, along with a small uncentering, can prevent 
the growth of high wavenumber noise. But in high resolution integrations (e.g., l°x 1.25° in 
the horizontal with 46 vertical layers), computational noise is present even with increased 
uncentering, especially at high latitudes and near orography (Moorthi, 1994). To alleviate 
this problem, the present implicit formulation of fourth order horizontal diffusion for a reg- 
ular latitude/longitude grid on the sphere has been developed. Some justification for the 
use of fourth order rather than second order diffusion is also presented. 

A detailed description of the solver is presented in Section 2. The justification for us- 
ing fourth order diffusion is discussed in Section 3. Results of applying the solver to a 
global shallow water semi-Lagrangian model are given in Section 4. Concluding remarks 
are presented in Section 5. The Fortran program of the solver is attached to the note as an 
Appendix. 


2 Description of The Solver 

2.1 Equations 


The fourth order horizontal diffusion equation takes the form 


and can be discretized as 


d<p 

dt 


-I<VU 


(f> n + l - 4> n 

At 


-tfv 4 ((i-o<r +, +e<r) 


(1) 

(2) 
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where £ is a parameter whose value is determined by the type of temporal discretization 
used; it is semi-implicit if f is ^ and implicit if £ is 0. (This terminology is adopted for 
consistency with the meteorological literature; the terms Crank-Nicholson and Euler implicit 
are corresponding descriptions used elsewhere.) 


Integrating (2) involves the solution of the following equation, in a manner similar to that 
used by Lindzen and Kuo (1969): 


vV +1 + 
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KAt (1 - 0 


<t> 


n + l 
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KAt (1 -0 


<t> n 


( 1-0 


vV 


( 3 ) 


i. e., 

V 4 <(>+r<t> = f (4) 

where r is a prescribed parameter, / is the forcing term and the superscript ( n + 1) has been 
omitted. Numerical solution of (4) on the sphere requires a conservative discretization of 
the V 4 operator, with a consistent treatment of the polar boundary conditions. Difficulties 
associated with these as well as with the cross derivative terms can be avoided by writing 
the equation in the form of two second order equations: 


where 


V 2 Z + r<j) = f 
V 2 <j>-Z = 0 


V 2 


1 d 2 1 8 ( 8 _\ 
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(6) 


2.2 Longitudinal and Latitudinal Discretization 


Since the coefficients in the above equations are independent of longitude, we can perform 
Fourier transforms using the Fast Fourier Transform (FFT) algorithm in the A— direction: 




( 7 ) 


Then the partial differential equations (5) and (6) lead to the ordinary differential equations 
(ODE) in vector form: 

LX+QX = F (8) 

where 
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The discretized forms of L and q at the interior points are: 

{Ly) J = ujy j+ 1 + Vjyj + 


sin 


9j = 
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for j = 1, 2, 3, . . . , N - 1 ( j = 0 and j — N represent the south and the north poles, 
respectively). We denote 


Aj = Wjl Bj — Vjl + Qj Cj — Ujl 


( 10 ) 


where I is the unit diagonal matrix. Thus, the discretized equation for the interior points, 
j = 1 , 2 , 3 , . . . , TV — 1 , can be written as: 


Ajjtj-x + BjXj + CjXj+i = Fj. 


( 11 ) 


2.3 Boundary Conditions at The Poles 

The boundary conditions at the south and the north poles are defined as: 

AqXq + Bo^i = (12) 

AwXn-i + B^Xn — Fjy. (13) 


At each pole, there should be a single unique value for the solution Z and <f > . For all 
longitudinal asymmetric components (k > 0), the boundary conditions for the Fourier 
components X at the poles simply reduce to: 

X 0 = 0 

Xn = 0 . 


The Fourier component ( k = 0) is the zonal mean, whose solution is derived from the ODE 
for it. The boundary conditions at the poles can be determined simply by discretizing (5) 
and (6) at the poles. We arrive at 
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for the south pole, and similarly 



for the north pole. The overbar denotes zonal average. These correspond to the boundary 
conditions of the form in (12) and (13) with the following coefficient matrices: 



2,4 Numerical Solution 

We again adopt the method of Lindzen and Kuo (1969) to solve the system (11) coupled 
with (12) and (13). Assume that the solution of (11) takes the following form: 

X, = aj x j+1 + a, ( 22 ) 
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where a 3 is a 2 x 2 matrix and (3j is a vector. We know from (12) that 


p 
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1 

& 

(24) 

If we substitute 

Xj-\ — Qtj— i X j + i 


into (11), we have 


X 3 — — + Bj) 1 CjXj + 1 + (Ajaj - 1 + Bj ) 1 ( Fj — 

AjPj-i) • (25) 

Comparing (25) and (22) leads to: 


a j = — (AjOtj-i + Bj) Cj 

(26) 

ft = (Ajaj-i + B,)~' [Pj - AjPj-i) 

(27) 

for j = 1, 2, 3, . . N — 1. With the following two equations: 


A^Xn-i 4- B^Xn = F/v 

(28) 

Xs-i = ckn-iXn + 

(29) 

we can obtain the solution for the north pole as: 


Xn - (Anoin-i + Bn) 1 (^Fn - AnPn- i) • 

(30) 

Now we recursively use (22) to obtain the complete set of Fourier components for each 
j = N - 1, N — 2, . . . , 2, 1, 0. Finally, we use inverse FFT to obtain the solution at grid 
points. 


2*5 Coding Strategy 

The coding strategy used here is similar to that of Moorthi and Higgins (1992, 1993). The 
FFT has been made very efficient by using the CRAY subroutine RFFTMLT to transform 
(both forward and backward) all latitudes simultaneously. In addition, the code for solving 
the second order O.D.E’s for the FFT components and the zonal mean is vectorized over the 
wave numbers. After vectorization, a speed-up by a factor of 40 was observed on the CRAY 
C-90 computer. The Fortran program with comprehensive documentation is supplied in the 
Appendix. 
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2.6 Extension to Higher Orders 


Extending the solver presented above to higher order is straightforward, both in formulation 
and in coding. Suppose we are to solve a m th order equation (m being an even number). In 
the place of (2.5) and (2.6), we now have a set of m/2 second order equations. The FFT, 
latitudinal discretization and numerical procedures are almost the same. The matrix size 
of 2x2 then becomes m/2 x m/2. As long as m is not too large, say, not larger than 8, the 
computational cost is not expected to increase dramatically. 


3 Justification for Using Fourth Order Diffusion 


It is well known that higher order diffusion can effectively damp the high wavenumber noise 
but leave the low wavenumbers very little affected. Here we present a discussion of the 
damping properties of the second and fourth order diffusion equations, in the context of 
both implicit and semi-implicit formulations. Though these properties have been discussed 
in earlier works (e.g., Mesinger and Arakawa, 1976), we present these here for the conve- 
nience of reference. We use the 1-D diffusion equation for illustration and then discuss the 
implication for spherical applications. 


The second order diffusion equation in its discretized form can be written as: 


<# +1 - <t>i 


At 


= m ((i - o (<w n+1 ) 7 + * ($r*<r)/) 


(31) 


where f = 0 corresponds to the implicit scheme and £ = | the semi-implicit scheme. The 
notation 6 XX is the finite difference approximation of the second order derivative, which is: 

A/+i + - 2.0A/ 

Vxx A — 2 

(Ax) 

Assume the solution takes the following form: 

= $ 0 AV fc/Ax . 

With second order centered difference for the spatial discretization, the damping rate, A 2 , 
becomes: 

1 — 4fr 2 sin 2 (— ) , v 

A ’ - 1 +4(l-fl A) ' <32) 

where r 2 = and m (= ttx) IS the wave number measured by grid intervals. 

The fourth order diffusion equation in its discretized form is 
1 


At — = -k ((i-o (^xx^h+a^-n/) 


(33) 
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where 


^ A + (^rx^)/_i 2.0 (S xx A)j 

Qxxxx™ — ~ To 

(Ax)' 

is the finite difference approximation of the fourth order operator. The damping rate, A 4l 


1 - 16£ r 4 sin 4 (^) 

l + 16 (l- 0 r 4 sin 4 (^) 


where r “ = (Sp* 


From (32) and (34) it is clear that for both second order and fourth order diffusion, the fully 
implicit treatment (f = 0) is unconditionally stable. When 0 < £ < 1, there is a critical 
value for K above which there is a reversal of the sign of the damping rate. When f 
there is also a critical value of K above which the scheme becomes unstable, this critical 
value being the smallest when the scheme is fully explicit (i.e., f = 1). When f is nonzero, 
it is possible to choose a value of K such that 2 Ax waves are completely eliminated. 

Figure 1 shows the damping rates as a function of wave number for the second and the fourth 
order diffusion equations in both implicit and semi-implicit formulations. The parameters 
r *2 and r 4 are so chosen that the 2Ax waves are completely eliminated with the semi- 
implicit formulations. Three features are observed: 1) In all cases, |A| < 1, which indicates 
unconditional stability; 2) The fourth-order equation is much more scale selective than the 
second order one; while the second order diffusion damps the 20Ax waves by about 5%, 
the fourth order leaves it almost unchanged; 3) The semi-implicit scheme can completely 
eliminate the 2Ax waves but the implicit scheme damps it only by a factor. 

The above results might mislead one to believe that the semi-implicit scheme is better than 
the implicit one since the former has the potential of eliminating the 2Ax waves completely 
with an appropriate diffusion coefficient. But for our case on the sphere with a uniform 
latitude/longitude grid, it may well be the situation that K is too large near the poles but 
too small near the equator. When K is very large, the 2Ax waves just change sign instead 
of being smoothed (see (34)). 

Figure 2a is the geopotential field of the initial conditon described in McDonald and Bates 
(1989), consisting of spherical harmonic wave number one. Figure 2b is the same field 
with 2Ax waves embedded. The result of applying the semi-implicit fourth order diffusion 
equation, with r 4 so chosen that the 2 Ax waves at the equator are eliminated, is shown in 
Figure 2c. The noise away from the equator remains. Careful examination shows that sign 
changes occur. Figure 2d represents the result after applying the implicit diffusion with 
the same value for the diffusion coefficient. Now the noise near the poles is damped, but 
some noise near the equator remains. Figure 2e is obtained by applying the semi-implicit 
formulation followed by the implicit formulation. It is almost identical to the original field 
in Figure 2a. Therefore, the use of the semi-implicit formulation of the diffusion equation 
alone should be avoided. The implicit formulation is preferred. If some noise near the 
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equator needs to be damped, the semi-implicit formulation, in addition to the implicit one, 
can be used if the additional computational cost is not a concern. 


4 Application to A Semi-Lagrangian Shallow Water Model 

In this section, we present some results of applying the fourth order diffusion to the shallow 
water model based on the momentum formulation, termed the (u,v) model as in Bates et 
al. (1994). This is the 2-D version of the NASA/GLA semi-Lagrangian GCM (Bates et a/., 
1993). 

4.1 Initial Condition of McDonald and Bates (1989) 

For a resolution of 512x257 on a uniform latitude/longitude grid and a timestep of one hour, 
a minimum value of 0.03 for the uncentering parameter is needed to suppress noise in the 
case of the geostrophically balanced initial condition of McDonald and Bates (1989). The 
initial geopotential is defined as: 

<f> (A, 9) = <f> + 2f2auosin 3 0cos 0sin A 

where fl and a are, respectively, the earth’s rotation and radius. The values <f> = 5.7879 X 10 
m 2 /s and u 0 = 20 m/s are used. The geopotential at 50 days of the integration is presented 
in Figure 3a. 

Now we replace the uncentering by the fourth order diffusion. In order to remove the residual 
noise near the equator, we apply the semi-implicit formulation in addition to the implicit 
one as discussed in the previous section. The diffusion coefficient for the semi-implicit 
formulation takes a value of 1.32 x 10 15 m 4 /s, which is just large enough to completely 
eliminate the 2Ax waves at the equator. The value needed for the implicit formulation is 
found to be dependent on the way the diffusion is applied. If the diffusion is applied to the 
geopotential only, a minimum value of 9.25 X 10 15 m 4 js is needed. The geopotential at 50 
days of the integration thus obtained is shown in Figure 3b. If the diffusion is applied only 
to the divergence, the minimum value needed is 11.825 X 10 15 m 4 /s. The geopotential at 50 
days of the corresponding integration is shown in Figure 3c. Remember that the divergence 
is not a basic prognostic model variable, rather it is an intermediate variable. This may be 
the reason why a relatively larger minimum value is needed. If we apply the diffusion to 
both the divergence and the geopotential, the minimum value for the diffusion coefficient 
becomes 5.28 X 10 15 m 4 /s , which is smaller than the value used in either of the previous 
two cases. The geopotential at 50 days of the integration is now presented in Figure 3d. 

The results shown in Figures 3b-d, which were obtained by applying the diffusion in various 
ways, demonstrate less damping of the long waves than is the case when using uncentering 
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(Figure 3a). Of the three cases (diffusion applied to geopotential alone, to divergence 
alone, and to both geopotential and divergence), it seems that the last two cases have the 
least undesirable damping effects. In addition, the results in Figures 3b-d all demonstrate 
the same phase speed, but the predicted geopotential of the uncentering case shows faster 
westward propagation of the long wave. It has been observed in our numerical results 
(not shown) that the westward phase speed monotonically increases as the uncentering 
parameter increases. This clearly demonstrates the negative impact of the uncentering on 
the long waves in both amplitude and phase. Improvement is obtained by replacing it by 
the fourth order diffusion. 

Examination of the evolution of the global mean geopotential indicates that the fourth order 
diffusion, regardless of the magnitude of the diffusion coefficient, has almost no impact on 
global mass conservation. 


4.2 Real Data Case 

The ECMWF analysis of 0000 UTC, 1 December 1992, is chosen. Again, the horizontal 
resolution is 512x257 points on a uniform latitude/longitude grid, but the timestep is chosen 
to be 30 minutes. The data are initialized using a digital filter (Lynch and Huang, 1992), 
with a cutoff period of 6 hours. The initialized geopotential is shown in Figure 4a. 

Without any diffusion, a value of 0.05 is needed for the uncentering parameter. The geopo- 
tential at day 3 of the integration is shown in Figure 4b. When we apply the fourth order 
diffusion to both geopotential and divergence, the model can be integrated stably without 
uncentering. The diffusion coefficient for the semi-implicit formulation is 1.32 x 10 15 m 4 /s, 
and that for the implicit formulation is 3.96 X 10 15 m 4 /s. However, there is still some resid- 
ual noise near the equator, which can be eliminated most effectively by using a small value 
of uncentering parameter (0.015 in this case). Figure 4c is the geopotential at day 3 of the 
integration thus obtained. The improvement of the representation of the resolved features 
in Figure 4c over Figure 4b is recognizable (The high in the upper left part of the model 
domain is stronger in Figure 4c). The advantages of using the fourth order diffusion should 
become more obvious when orography is present, in which case a much larger value for the 
uncentering parameter (undesirable) is otherwise needed to control the high wavenumber 
noise in mountainous regions. 


5 Application to a 3-D Semi-Lagrangian GCM 

We have incorporated the fourth-order horizontal diffusion into the NASA/GLA semi- 
Lagrangian semi-implicit (SLSI) general circulation model(GCM) with full physics. The 
adiabatic part of this GCM is described in Bates et al. (1993) and the version with full 
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physics is described in Moorthi et al (1994). Here, we present some results to illustrate 
the use of the fourth-order diffusion. 

The version of the model used here has twenty sigma layers in the vertical and a horizontal 
resolution of 2 degree latitude by 2.5 degree longitude. The initial state for the experiments 
is taken from the NASA/GLA Data Assimilation Office analyses for November 1, 1992, 
performed during the 1992 fall SPADE mission, using the GLA Eulerian GCM. The fourth 
order diffusion in its implicit formulation is applied to both temperature and moisture. The 
diffusion coefficient is taken a s K = min ^ Q - 25 ^ 1Ql5 , 5.0 X 10 15 ). Here cq is the value of the 

vertical coordinate sigma at the middle of the I th layer. This choice is made particularly 
due to the fact that there exists stronger undesired noise in the model integration above the 
100 hPa level than at lower levels. The associated e-folding time for the shortest resolvable 
waves at the equator is approximately 2 days near the bottom and 0.1 day near the top. 

We have performed four ten-day runs, two without diffusion but with values of 0.05 and 
0.07 for the uncentering parameter, and two with diffusion and with a value of 0.05 for the 
uncentering parameter. Both runs without diffusion show significant noise at the upper 
level of the model. An example is shown in Figure 5, where the temperature field at 50 
hPa level is shown for the forecasts at days 3, 5, 7, and 10 respectively for the run with 
the uncentering parameter value of 0.05. It can be seen that the temperature field is noisy 
and the noise progressively increases near the equator as the integration proceeds. The 
corresponding field for the forecast with diffusion is shown in Figure 6. It can be seen that 
the inclusion of diffusion has effectively reduced the high wavenumber noise. Of course, one 
can obtain even better smoothing by increasing the value of the diffusion coefficient. The 
results for the uncentering parameter value 0.07 are similar to those shown above. Because 
the noise in the upper level equatorial region (which may be the result of the gravity waves 
generated by convective processes), cannot be efficiently eliminated through judicial choice 
of the uncentering parameter, inclusion of a higher order diffusion becomes very important. 


6 Concluding Remarks 


A direct solver for implicit fourth order horizontal diffusion on the sphere has been devel- 
oped. It has been made efficient by vectorizing over the FFT wave spectrum. Extension 
of the formulation to the solution of even higher order diffusion is straightforward, and the 
additional coding effort needed is trivial. 

Applying the solver to a global shallow water semi-Lagrangian model leads to positive 
results. It gives a more accurate representation of larger scales in both amplitude and 
phase than does the use of the uncentering method of damping. It is found that horizontal 
diffusion has almost no impact on the global mass conservation. The solver is also applied 
to the 3-D NASA/GLA semi-Lagrangian semi-implicit GCM and found to be helpful in 
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controlling the high wavenumber noise. The solver takes approximately 4% of the total CPU 
time on a single processor of the CRAY C98 computer when applied to the temperature 
field in the adiabatic version of the 3-D GCM. 

Acknowledgements. This research was supported by the NASA Global Atmospheric Mod- 
eling and Analysis Program under Grant 578-41-16-20 
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APPENDIX A 

Fortran Program of The Solver 


Pag* 1 


SUBROUTINE DIRS0L4(IMA, IM, JM,AE, INIT1, INIT2 , DIFF, FI, FO, XK, DT, IMPL) 

C Subroutine DIRS0L4 

Q★**★★★*****★***★***★★*■**•****★^r**^^***★★*★*** , *■*★★★*** + + ** + * + * , + ***■*'*** , *‘*■**'*■* 

C WRITTEN BY: Yong Li, GSFC Code 910. 3 , Apr 1994 

C 

C PURPOSE: This routine solves the fourth-order diffusion equation on the 

C sphere implicitly or semi-implicitly (Crank-Nicolson scheme), that is, 

C to solve the following equation: 

C 

C (\partial phi over \partial t) = - K \delta*4 (phi) 

C 

C where K is the diffusion coefficient. 

C 

C Implicit Soluton: (phi* (n+1 ) -phi* (n) ) /DT = - K \del*4 (phi* (n+1)) 

C 

C Semi- Implicit : (phi* (n+ 1 ) -phi* (n ) ) /DT = 

C - K (1/2) \del*4 (phi* (n+1) +phi* (n) ) 

C 

C Each of the above discrete equation converts to the solution 

C of a linear fourth-order partial differential equation on the shpere 
C of the following form: 

C 

C \delta*4 (phi)* (n+1) - r*phi*(n+l) = f 

C 

C The strategy of solution is similar to that in Lindzen and Kuo (1969, 

C Mon. Wea. Rev., p732-734) , which is summarized as: 

C 

C 1) . Convert the \del*4 -equation into a set of two \del*2 equations. 

C Now the model variables become a vector of 

C ( \del*2 (phi ) * (n+1) , (phi)* (n+1) ) ; 

C 

C 2) . Since the coefficients are only latitude-dependent, FFT along the 

C latitudes are performed to transform the second-order vector PDE 

C to a series of second-order vector ODE'S for the FFT coefficients; 

C 

C 3). Solve the set of second-order ODE's; 

C 

C 4) . Perform inverse FFT to obtain the solution (phi)* (n+1) on gridpoints. 

C 

C 

C CALLED FROM: Any routine which applies this smoother. 

C 

C SYSTEM ROUTINES USED: 

C FFTFAX : Cray routine which initializes the FFT routines; 

C RFFTMLT: Cray FFT routine (Cray Manual: SR-2081 8.0) 

C 

C SUBROUTINES CALLED: 

C CAL_R: calculate the coefficient r from K and and DT; 

C FORCE: calculate the forcing term f for the \del*4 equation; 

C S20DE: solve the set of the second-order ODE's for FFT 

C coefficients k > 0; 

C SOMEAN: solve for the zonal mean; 

C WRAPD: wrap ghost points along the latitudes. 

C 

C VARIABLES: 

C 

C INPUT : 

C AE-radius of earth (meters) ; 

C DIFF-logical variable meaning discrete FFT when it is .TRUE.; 

C DT-time step (seconds) 

C FI-input field before the smoothing (phi)*(n); 

C IM-number of grid points along latitudes; 

C IMA-dimension in X-direction (=IM + # of ghost points); 

C IMPL -fully implicit scheme ( = 1) and semi-implicit scheme (=0) ; 

C INITl-logical variable (need to be TRUE) ; 

C INIT2-logical variable (need to be TRUE) ; 

C XK-diffusion coefficient K (m*4/s) . 
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OUPUT : 

FO-output field after smoothing (phi)~(n+l). 


LOCAL: 

A-array storing the FFT coef f icinets ; 

AM-coefficient array for the 2nd-order ODE (for the zonal mean) 
at j-1 point; 

AMK-coef f icient array for the 2nd-order ODE for the FFT coefficients 
at j-1 point; 

B-working array for the FFT; 

BM-coeff icient array for the 2nd-order ODE (for the zonal mean) 
at j point; 

BMK-coeff icient array for the 2nd-order ODE for the FFT coefficients 
at j point; 

CM-coeff icient array for the 2nd-order ODE (for the zonal mean) 
at j + 1 point; 

CMK-coeff icient array for the 2nd-order ODE for the FFT coefficients 
at j+1 point; 

CPH-COS ( theta) where theta is the latitude; 

DNl-real part of the forcing term for the 2nd-order ODE or forcing 
for the equation for the zonal mean; 

DNIK-real part of the forcing term for the 2nd-order ODE or forcing 
to the equation for the FFT coefficients; 

DN2K-imaginery part of the forcing term for the 2nd-order ODE or forcing 
to the equation for the FFT coefficients; 

F-array storing the forcing term for the \del A 4 equation; 

FM-zonal mean of the input field (phi) A (n); 

FNl-the solution to the 2nd-order ODE for the zonal mean; 

FNIK-real part of the solution to the 2nd-order ODE for the FFT 
coefficients ; 

FN2K-imaginery part of the solution to the 2nd-order ODE for the FFT 
coefficients ; 

IX-working array for the FFT routine; 

TR-working array for the FFT routine. 

U-coef f cients for intermediate calculation; 

V-coef f cients for intermediate calculation; 

W-coeff cients for intermediate calculation; 

WVSQ-array storing coefficients for the X-direction derivative term 
in the \del A 2 equation; 


INTEGER FORWARD, BACKWARD 
PARAMETER (FORWARD = -1, BACKWARD = 1) 

DIMENSION FI ( IMA, JM) , FO(IMA,JM) 

DIMENSION IX(IM), TR ( 3 * ( IM+2 ) + 1 ) , A ( IM+2 , JM-2 ) , B(IM*2,JM-2) 


DIMENSION AM(2, 2, JM) , BM(2,2,JM), CM(2,2,JM), F(IMA,JM) 

#, DN1 { 2 , JM) , DN2 ( 2 , JM) , FN1(2,JM), FN2(2,JM) 

#, FM(JM), CPH ( JM, 2 ) , U ( JM) , V(JM) , W(JM) 

#, WV { IM/ 2 ) , WVSQ ( IM/2 ) 

DIMENSION AMK ( IM/2 , 2,2, JM) , BMK ( IM/2 , 2 , 2 , JM) , CMK ( IM/2 , 2 , 2 , JM) 
#, DN1K { IM/2 , 2 , JM) , DN2K ( IM/2 , 2 f JM) 

#, FN1K { IM/2 , 2 , JM) , FN2K ( IM/ 2 , 2 , JM) 


LOGICAL INIT1, INIT2 , DIFF, FIRST 

PI = ACOS (-1.0) 

TWOPI = 2.0 * PI 
PI02 = PI/2.0 
IMB2 = IM/2 
I JM = (IMB2-1) 

I JMM = I JM*22 + JM 

IJM2 = I JMM + (IM+2) * (JM-2) 
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CALL CAL_R ( XK , DT , R , IMPL ) 

CALL FORCE { IMA, IM r JM, FI , F, AE , R, IMPL) 

CALL WRAPD ( IMA , IM , JM , F ) 

IF (INIT1) THEN 

CALL FFTFAX ( IM, IX, TR) 

JMM1 = JM - 1 

JMM2 = JM - 2 

IMP1 = IM + 1 

IMP2 = IM + 2 

LEN = IMB2 - 1 

FIM = 1.0 / FLOAT (IM) 

DLM = TWOPI *FIM 
DPH = PI/ FLOAT (JMM1) 

RDLM = 1.0 /DLM 
RDPH = 1.0 /DPH 

DY = AE*DPH 
AESQ = AE*AE 

DO J = 1, JM 

TEM = - PI02 + { J-l ) *DPH 
CPH { J , 1 ) = COS (TEM+0 . 5*DPH) 

CPH ( J , 2 ) = COS (TEM) 

ENDDO 

CPH (1,2) =0.0 
CPH ( JM, 2 ) = 0.0 

FIRST = .TRUE. 

INIT2 = .TRUE. v 

ENDIF 

IF (INIT2) THEN 
DO J = 2, JMM1 

U(J) = CPH ( J , 1) /CPH( J, 2) 

W<J) = CPH ( J- 1 , 1 ) /CPH ( J , 2 ) 

V(J) = - 1.0 * (U(J)+W(J) ) 

ENDDO 

IF (DIFF) THEN 
DO I = 1, I MB 2 

WVSQ(I) = (SIN ( 0 . 5 * I * DLM ) * (RDLM* 2 . 0) ) **2 
ENDDO 
ELSE 

DO I = 1, IMB2 
WVSQ(I) = 1*1 
ENDDO 
ENDIF 

ENDIF 

A ( : , : ) =0.0 

DO J = 1, JM 
FM( J) = 0.0 
DO I = 1, IM 

FM ( J ) = FM ( J ) + F ( I , J ) 

ENDDO 

ENDDO 

DO J = 1, JM 

FM ( J) = FM ( J) *FIM 
ENDDO 

DO J = 2, JMM1 
DO I = 1, IM 

A ( I , J-l ) = F ( I , J) - FM ( J) 

ENDDO 
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A ( IM+ 1 , J - 1 ) = A ( 1 , J-l ) 

A { IM+2 , J- 1 ) = A ( 2 , J-l ) 

ENDDO 

CALL RFFTMLT { A , B , TR , IX, 1, IMP2 , IM, JMM2 , FORWARD) 

A ( 1 , : ) =0.0 
A ( 2 , : ) = 0.0 

A ( IMP2 , : ) =0.0 


DO J = 2 , JMM1 
DO I = 1, IMB2 

AMK ( I , 1 , 1 , J ) = W{J) 

AMK ( I , 2 , 2 , J ) = W(J) 

AMK ( 1 , 1 , 2 , J) =0.0 
AMK<I,2,1,J) =0.0 

CMK { I , 1 , 1 , J) = U ( J ) 

CMK < I , 2, 2, J) = U{J) 

CMK { I , 1,2, J) = 0.0 
CMK ( I , 2 , 1 , J ) = 0.0 

BMK ( 1 , 1 , 1 , J ) = - <DY**2) *WVSQ{I) / <AE*CPH(J, 2) ) **2 

# + V(J) 

BMK ( I , 1 , 2 , J ) = R* (DY**2 ) 

BMK(I, 2, 1, J) = - DY* *2 

BMK { 1 , 2 , 2 , J ) = - (DY* * 2 ) *WVSQ ( I ) / { AE*CPH { J , 2 ) ) **2 

# + V(J) 

DN1K ( I , 1 , J ) = (DY**2) *A(I+I+1, J-l) 

DN1K ( I , 2 , J ) =0.0 

DN2K ( I , 1 , J ) = <DY**2) *A{I+I+2 # J-l) 

DN2K ( I , 2 , J ) =0.0 

ENDDO 
ENDDO 


DN1K ( : , : , 1) = 0.0 
DN2K = 0.0 


AMK ( : 

,1,1,1) 

= 

1.0 

AMK { : 

,2,2,1) 

= 

1.0 

AMK { : 

,1,2,1) 

= 

0 . 0 

AMK ( ; 

,2, 1, 1) 

= 

0.0 

BMK { : 

, : , : , 1) 

= 

0 . 0 

DN1K ( 

: , : , JM) 

- 

0.0 

DN2K { 

: , : , JM) 

= 

0.0 


AMK ( JM) = 0.0 


BMK { : f 1, 1, JM) = 1.0 
BMK { : ,2,2, JM) = 1.0 
BMK ( : , 1, 2, JM) = 0.0 
BMK ( : ,2, 1, JM) = 0.0 

CALL S2 ODE { IMB2 , DN1K, DN2K , FN IK , FN2K , AMK , BMK , CMK , JM ) 

DO J = 2, JMM1 
DO I = 1, IMB2 

Ad+I+l, J-l) = FN1K ( 1 , 2 , J ) 

A { 1 + 1 + 2 , J-l) = FN2K ( 1 , 2 , J) 

ENDDO 

ENDDO 


DO J = 2, JMM1 
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AM ( 1 , 1 , J) 
AM { 2 , 2 , J) 
AM<1, 2, J) 
AM( 2 , 1 , J) 

CM ( 1 f 1 , J ) 
CM { 2 , 2 , J ) 
CM{1, 2, J) 
CM { 2 , 1 , J ) 

BM { 1 , 1, J) 
BM < 1 , 2 , J ) 
BM { 2 , 1, J) 
BM ( 2 , 2 , J ) 

DN1 { 1 , J) 
DN1 (2, J) 
ENDDO 


= W(J) 

= W(J) 

= 0.0 
= 0.0 

= U(J) 

= U(J) 

= 0.0 
= 0.0 

= + V ( J ) 

= R* ( DY* *2 ) 

= - DY**2 

= BM { 1 , 1 , J ) 

= (DY**2 ) *FM ( J) 

= 0.0 


DN1 (1,1) = <DY**2)*FM<1) 
DN1 (2,1) = 0.0 


AM{1, 1, 1) 
AM(1, 2, 1) 
AM (2,1,1) 
AM {2,2,1} 


- 4.0 

R* (DY* *2 ) 

- (DY* *2 ) 

- 4.0 


BM (1,1,1) = 4.0 
BM (1,2,1) = 0.0 
BM (2,1,1) = 0.0 
BM (2,2,1) = 4.0 


DN1 ( 1 , JM) = (DY**2) *FM(JM) 

DN1 { 2 , JM) = 0.0 

AM {1,1, JM) = 4.0 
AM (1,2, JM) = 0.0 
AM (2,1, JM) = 0.0 
AM (2,2, JM) = 4.0 

BM {1,1, JM ) = - 4.0 
BM (1,2, JM) = R* ( DY* *2 ) 

BM {2,1, JM ) = - (DY**2) 

BM (2,2, JM) = - 4.0 

FN1 { : , : ) = 0.0 

CALL SOMEAN (DN1 , FNl , AM, BM, CM, JM) 

CALL RFFTMLT { A , B , TR , IX, 1, IMP2 , IM , JMM2 , BACKWARD ) 

DO J = 1, JMM2 
DO I = 1, IM 

FO ( I , J+l ) = A ( I , J ) + FN1 ( 2 , J+ 1 ) 

ENDDO 

ENDDO 

DO 1 = 1, IM 

FO(I,l) = FNl (2,1) 

FO ( I , JM) = FNl (2 , JM) 

ENDDO 

CALL WRAP1 ( IMA, IM, JM, FO ) 

RETURN 

END 

SUBROUTINE CAL_R ( XK , DT , R , IMPL ) 

AB = ( 1 . 0+FLOAT < IMPL) ) *XK*DT 
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R = 2 .O/AB 

RETURN 

END 

SUBROUTINE DEL2 ( IMA, IM, JM, A, DA, AE) 

C Subroutine DEL2 

C****************************************************************** 

C WRITTEN BY: Yong Li f GSFC code 910.3, Apr 1994 
C 

C PURPOSE: calculate the \del A 2(A) on the sphere. 

C 

C CALLED BY: 

C DEL4 

C 

C SUBROUTINES CALLED: None 

C 

C VARIABLES: 

C 

C INPUT : 

C A- input field; 

C AE-radius of earth; 

C IM-number of total grid points in X-direction; 

C IMA -dimension size in X-direction; 

C JM-number of total grid points in Y-direction; 

C 

C OUTPUT : 

C DA-\del*2 of A 

dimension a(ima,jm), da(ima,jm), dx(jm) 
dimension al (500) , a3{500), bp<500), bm{500) 

impl = im + 1 

imp2 = im + 2 

pi = acos (-1.0) 

twopi = 2 . 0*pi 

dim = twopi/f loat ( im) 

dph = pi/float (jm-1) 

dy = ae*dph 

a2=l . 0d0/dy**2 

do j = 2, jm-1 

theta = - pi*0.5 + f loat ( j -1 ) *dph 
bp(j) = cos (theta+0 . 5 *dph ) /cos ( theta ) 
bm( j ) = cos ( theta-0 . 5*dph) /cos ( theta) 
dx(j) = ae*dlm*cos ( theta) 
al< j) = 1 . 0/dx ( j ) **2 
a3(j) = 2 . 0*al { j ) +a2* (bp ( j ) +bm( j ) ) 
enddo 


do 20 i=2, impl 
do 30 j =2 , jm-1 

da ( i, j ) =al ( j ) * (a (i+1, j ) +a (i-1, j ) ) 

# + a2 * ( bp ( j ) * a ( i , j + 1 ) +bm( j ) *a(i, j -1) ) 

# -a3 ( j ) *a ( i, j ) 

30 continue 

20 continue 


C 

C 

C 


do j = 2, jm-1 

da ( 1 , j ) = da ( im+ 1 , j ) 
do i = im+2, ima 

da ( i , j ) = da ( i - im , j ) 
enddo 
enddo 


Laplace at the south pole 


AAC =0.0 
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do 40 i = 2, impl 

AAC = AAC + a ( i , 2 ) 

40 continue 

da(l, 1) = { 2 . 0 /dy) ** 2 * (AAC/ float { im) -a ( 1 , 1 ) ) 

do 50 i = 2 r ima 
da ( i , 1 ) = da ( 1 , 1 ) 

50 continue 


Laplace at the north pole 

AAC =0.0 

do 60 i = 2 f impl 

AAC = AAC + a(i, jm-1) 

60 continue 

da ( ima, jm) = (2 . 0/dy) ** 2 * (AAC/ float ( im) -a ( 1, jm) ) 

do i = 1, ima- 1 

da(i,jm) = da(ima # jm) 
enddo 

return 

end 

SUBROUTINE DEL4 ( IMA, IM, JM, A, DA, AE) 

*★★★■**★**•*****★★★****■**★**★*******★**•#*★★★*★★*************★**★***★*★**■*■* 
Subroutine DEL4 

****★**★**★*★*★★★*★**★*******★*★*★***★*★★***★*★***★★*★**★*★*★**★★*****★* 

WRITTEN BY: Yong Li, GSFC code 910.3, Apr 1994 

PUEPOSE: calculate the \del /v 4(A) by applying a conservative \del JN 2 twice. 
CALLED FROM: 

FORCE: calculate the forcing term 

SUBROUTINES CALLED: 

DEL2 : \del A 2 operator; 

WRAPD: wrap up the ghost points in X-direction. 

VARIABLES : 

INPUT: 

A-input field; 

AE-radius of earth; 

IM-number of total grid points in X-direction; 

IMA-dimension size in X-direction; 

JM-number of total grid points in Y-direction. 

OUTPUT : 

DA-\del"4 (A) 

****»*****************i.*****lr-»**********************i******** ************** 

DIMENSION A (IMA, JM) , DA ( IMA, JM) , B(IMA,JM) 

CALL DEL2 ( IMA, IM, JM, A, B, AE) 

CALL WRAPD( IMA, IM, JM, B) 

CALL DEL2 ( IMA, IM, JM, B, DA, AE) 

CALL WRAPD ( IMA, IM,JM, DA) 

RETURN 

END 

SUBROUTINE FORCE (IMA, IM, JM, A, F, AE, R, IMPL) 

**************. **************. w********************************* 

Subroutine FORCE 

*******************************************************«r*Hr****** 

WRITTEN BY: Yong Li, GSFC code 910.3, Apr 1994 
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PURPOSE: calculate the forcing term for the \del~4 -equation . 
CALLED FROM: DIRS0L4 
SUBROUTINES CALLED: 

DEL4 : calculate the \del A 4 operator of field A; 

WRAPD: wrap up the ghost points in X-direction; 

VARIABLES : 


INPUT: 

A-array storing the field to be smoothed; 

AE-radius of earth; 

IM-number of total grid points in X-direction; 

IMA-dimension size in X-direction; 

IMPL-for implicit (=1) , for semi-implicit (=0) ; 

JM-number of total grid points in Y-direction; 

R-coef f icient for the zero-order term in the \del A 4 equation; 

OUTPUT : 

DA-array storing \del A 4(A); 

F-array storing the forcing term for the \del A 4 equation 
****** + *** + **************** + ******************* + **** + *** + * * *#***★***★*★*★ 

DIMENSION A {IMA, JM) , DA ( IMA , JM ) f F(IMA,JM) 

IF ( IMPL . EQ . 1 ) THEN 
F ( : , : ) = R*A(:, :) 

ELSE 

CALL DEL4 { IMA , IM , JM , A , DA , AE ) 

CALL WRAPD ( IMA , IM , JM r DA ) 

F ( : , : ) = R*A ( : , :) - DA ( : , : ) 

ENDIF 


RETURN 

END 


SUBROUTINE INVERSE 1 (A, B) 


Subroutine INVERSE1 


WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 

PURPOSE: Find the inverse of a 2x2 matrix 

VARIABLES : 

INPUT: 

A-input matrix to be inverted; 

OUTPUT : 

B-inverse of matrix A; 

DIMENSION A (2, 2) , B(2,2) 

DD = A ( 1 , 1) *A ( 2 , 2 ) -A ( 1 , 2) *A ( 2 , 1) 

B ( 1 , 1 ) = A ( 2 , 2 ) /DD 

B(l, 2) = - A(l,2) /DD 
B ( 2 , 1 ) = - A(2, 1) /DD 
B ( 2 , 2 ) = A(l, 1) /DD 

RETURN 
END 

SUBROUTINE INVERSE2 (LEN, A, B) 

C Subroutine INVERSE2 

Q*»«***««t****«*t*t*tt****t*t*tt**ttt««***t*****»tt1r*t******t**** 

c 

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 
C 
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C PURPOSE: Find the inverse of a a number (LEN) of 2x2 matrix 
C 

C VARIABLES: 

C INPUT : 

C A- input of LEN 2x2 matrices to be inverted simultaneously; 

C OUTPUT : 

C B-inverse the LEN 2x2 matrices A; 


Dimension A(LEN,2,2), B(LEN,2,2), DD(LEN) 

DO KK = 1, LEN 

DD<KK) = A(KK, 1, 1) *A(KK,2, 2) -A(KK, 1,2) *A(KK, 2, 1) 
ENDDO 


DO KK = 1, LEN 
B { KK ,1,1) = 
BUCK, 1,2) = 
BUCK, 2,1) = 
BUCK, 2,2) = 
ENDDO 


A ( KK ,2,2)/ DD { KK ) 

- A (KK, 1,2) /DD (KK) 

- AUCK, 2, 1) / DD ( KK ) 
A (KK, 1,1) /DD (KK) 


RETURN 

END 

SUBROUTINE MULTAB ( A, B, C, N) 

***★**★★**★****★**★**★★★*★★*★★*★*★★★★★*****+*★★**★****★***★**★* 
Subroutine MULTAB 


WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 

PURPOSE: Multiplication of two NxN matrix C = A B 

VARIABLES : 

INPUT: 

A-input matrix; 

B-input matrix; 

OUTPUT : 

C-output matrix. 

DIMENSION A (N,N) , B (N, N) , C(N,N) 

DO I = 1, N 
DO J = 1, N 

C(I,J) = 0.0 
DO K = 1, N 

C ( I , J ) = C(I,J) + A(I,K> *B(K, J) 

ENDDO 

ENDDO 

ENDDO 

RETURN 

END 

SUBROUTINE MULTABK (LEN, A, B, C, N) 

£*★*****★****************★***★★*■*★***★★★**★*★***★*****★*★★★*★★★★ 
C Subroutine MULTABK 

£*****★***************★*•*★**★****★★★★***★★★★★★★★*★★★★*********** 

C 

C WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 
C 

C PURPOSE: Multiplication of a series of two NxN matrix C = A B 
C 

C VARIABLES: 

C INPUT : 

C A-input of LEN 2x2 matrices; 

C B-input of LEN 2x2 matrices; 

C OUTPUT : 

C C-output matrix. 

DIMENSION A(LEN,N,N) , B (LEN, N, N) , C (LEN,N, N) 

#, WRKl(LEN), WRK2 (LEN) 


21 



oooonoooooooooo noonnonnnononno 


Page 10 


DO I = 1, N 

DO J = 1, N 

C< : , I, J) = 0.0 
DO K = 1, N 

WRK1 ( : ) = C< : , I, J) 

DO KK = 1, LEN 

WRK2 (KK) = A (KK, I , K) *B (KK, K, J) 

ENDDO 

C ( : , I , J) = WRK1 { : ) + WRK2 ( : ) 

ENDDO 

ENDDO 

ENDDO 

RETURN 

END 

SUBROUTINE MULTAV { A, V, W, N) 

Subroutine MULTAV 

**#******** + ***** + ********* + **** + *-*-#*** + **-** + **************** + * 

WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 

PURPOSE: Multiplication of an NxN matrix A and a vector V: W = A V 

VARIABLES : 

INPUT: 

A-input matrix; 

B- input vector; 

OUTPUT: 

C-output vector. 


DIMENSION A(N,N), V(N), W(N) 

DO I = 1, N 
W ( I ) = 0.0 
DO J = 1, N 

W(I) = W(I) ♦ A ( I , J) *V(J) 

ENDDO 

ENDDO 

RETURN 

END 

SUBROUTINE MULTAVK ( LEN r A , V , W , N ) 

fr***********************************************************-*** 

Subroutine MULTAVK 

***** + *********** + ****-*********-**-*#**************-****■**** + * + **# 

WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 

PURPOSE: Multiplication of an NxN matrix A and a vector V: W = A V 

VARIABLES : 

INPUT: 

A-input of LEN 2x2 matrices; 

B- input of LEN vectors; 

OUTPUT : 

C-output vector. 

DIMENSION A(LEN,N,N) , V { LEN, N) , W (LEN, N) , WRK(LEN) 

DO I = 1, N 
W ( : , I ) =0.0 
DO J = 1, N 

WRK { : ) = W ( : , I ) 

W ( : , I ) = WRK(:) + A < : , I , J) *V < : , J) 

ENDDO 

ENDDO 

RETURN 

END 
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SUBROUTINE S20DE (LEN, DN1 , DN2, FN1 , FN2 , AM, BM, CM, JM) 
*******-****-***•**********■************* + ***-********** + ***«********** 

Subroutine S20DE 

****+**************+******+***********.*#************************** 

WRITTEN BY: Yong Li, GSFC Code 910.3, Apr 1994 

PURPOSE: This routine solves the second-order ODE which originally results from 
the \del A 4 equation on the sphere. But appplication is not so restricted 
as far as the discretized equation is in the format as (1) . 

The unknown to be solved is a 2x1 vector {Lindzen and Kuo, 1969, MWR) . 

Therefore this routine is used here to solve for the FFT coefficients (k > 0) 
in the \del A 4 direct solver if the the B.C's are chosen to be zero for the poles. 

The discretized form of the ODE is: 

A(j) X(j-l) + B ( j ) X(j) + C(j) X(j + 1) = F< j ) (1) 

(for j = 2, 3, 4, JM-1) 

The boundary conditions are: 

A ( 1 ) X(l) + B ( 1 ) X { 2 ) = F ( 1) (2.1) 

and A(JM) X(JM-l) + B(JM) X(JM) = F(JM) (2.2) 

The coefficient A, B and C are all real. The forcing F and hence the 
solution X have both real and imaginery parts. 

CALLED FROM: 

DIRSOL4 


SUBROUTINES CALLED: 

INVERSE1 : inversion of one 2x2 matrix; 

INVERSE2 : inversion of a number of 2x2 matrices; 


MULTAB: 
MULTABK : 
MULTAV : 
MULTAVK 


multiplication of two matrices; 
multiplication of a number of two matrix set; 

a matrix and a vector . 
a number of a matrix-vector sets . 


multiplication of 
multiplication of 


VARIABLES : 


INPUT: 

AM-array storing the coefficient matrix A as in (1); 

BM-array storing the coefficient matrix B as in (1) ; 

CM-array storing the coefficient matrix C as in (1); 

DNl-array storing the real part of the vector of the forcing 
term F in ( 1) ; 

DN2 -array storing the imaginery part of the vector of the forcing 
term F in (1) ; 

JM-number of total grid points including the two boundary points . 
OUTPUT: 

FNl-array storing the real part of the solution vector X of (1) ; 
FN2-array storing the imaginery part of the solution vector X of (1); 

LOCAL : 

ALPHA-array storing intermediate coefficient matrix alpha 
(Lindzen and Kuo, 1969); 

BETAl-real part of vector beta (Lindzen and Kuo, 1969); 

BETA2- imaginery part of vector beta (Lindzen and Kuo, 1969); 

VRK1 -working array; 

VRK2 -working array; 

VRK3 -working array; 

VRK4 -working array; 

VRK5 -working array; 

VRK6 -working array; 

WRK1 -working array; 

WRK2 -working array; 

WRK3 -working array; 

WRK4 -working array; 
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DIMENSION DN1( LEN, 2, JM) , DN2 (LEN, 2 , JM) , FN1 (LEN, 2 , JM) 

#, FN2 (LEN, 2, JM) 

DIMENSION AM(LEN, 2, 2, JM) , BM (LEN, 2 , 2 , JM) , CM (LEN, 2, 2, JM) 
DIMENSION ALPHA(LEN, 2,2, JM) , BETA1 (LEN, 2 , JM) , BETA2 (LEN, 2 , JM) 
DIMENSION WRK1 (LEN, 2, 2) , WRK2 (LEN, 2 , 2 ) , WRK3 (LEN, 2 , 2) 


#, 

WRK4 (LEN, 2, 2) 



# f 

VRK1 (LEN, 2) , VRK2 ( LEN , 2 ) , 

VRK3 (LEN, 2 ) , VRK4<LEN,2) 

#, 

VRK5 (LEN, 2) , VRK6(LEN,2) 


WRK1 { : , : , 

;) = AM( 1) 



CALL INVERSE2 (LEN, WRK1, WRK2 ) 



WRK1 ( : , : , 

: ) = - WRK2 (:,:,:) 



WRK4 ( : , : , 

:) = BM (:,:,:,!) 



CALL MULTABK (LEN, WRK1, WRK4 , 

WRK3 , 

2) 

ALPHA ( : , : 

, : r l) = WRK3<:, 



VRK1 { : , : ) 

= DN1{ : , : , 1) 



VRK2 < : , : ) 

= DN2 ( : , : , 1 ) 



CALL MULTAVK (LEN, WRK2 , VRK1, 

VRK3 , 

2) 

CALL MULTAVK (LEN, WRK2 , VRK2 , 

VRK4, 

2) 

BETA1 { : , : 

,1) = VRK3 ( : , ;) 



BETA2 ( : , : 

,1) = VRK4 ( : , :) 



DO J = 2, 

JM- 1 



WRK1 { : 

, : , : ) = ALPHA ( : , : , : , J 

-1) 


WRK2 ( : 

, :> = AM( J) 




CALL MULTABK (LEN, WRK2 , WRK1, WRK3 , 2) 

WRK1 (:,:,:) = WRK3 (:,:,:) +BM(:,:,:,J) 

CALL INVERSE2 ( LEN , WRK1 , WRK2 ) 

WRK1 ( : , :, :) = CM< J) 

CALL MULTABK (LEN, WRK2 , WRK1, WRK3 , 2) 

ALPHA ( J) = - WRK3 (:,:,:) 

WRK1 (:,:,:) = AM ( : , : , : , J) 

VRK1(:,:) = BETA1 ( : , : , J-l) 

VRK2 ( : , : ) = BETA2 ( : , : , J- 1 ) 

CALL MULTAVK ( LEN , WRK1, VRK1, VRK3 , 2) 
CALL MULTAVK ( LEN , WRK1, VRK2 , VRK4 , 2) 

VRK3 ( : , : ) =DN1(:,:,J) - VRK3 ( : , : ) 

VRK4 ( : , : ) =DN2(:,:,J) - VRK4 ( : , : ) 

CALL MULTAVK (LEN, WRK2, VRK3 , VRK1, 2) 
CALL MULTAVK (LEN, WRK2 , VRK4 , VRK2, 2) 

BETA1 ( : , : , J) = VRK1 ( : , : ) 

BETA2 ( : , : , J) = VRK2 ( : , : ) 

ENDDO 


VRK1 ( : , 

) = DNl ( : , : 

/ JM) 

VRK2 ( : , 

) = DN2 ( : , : 

, JM) 

VRK3 { : , 

) = BETA1 ( : 

, : , JM- 1 ) 

VRK4 { : f ; 

) = BETA2 ( : 

, : , JM- 1 ) 

WRK1 { : , ; 

:, :) = AM( : , 

i , JM) 
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CALL MULTAVK ( LEN , WRK1, VRK3 , VRK5 , 2) 
CALL MULTAVK < LEN , WRK1, VRK4 , VRK6 , 2) 

VRK1 ( : , : ) = DN1 ( : , : , JM) - VKK5 ( : , : ) 

VRK2 ( : , : ) = DN2 ( : , : , JM) - VRK6 ( : , : ) 

WRK1 = ALPHA (:,:,:, JM-1 ) 

WRK2 (:,:,:) = AM JM) 

CALL MULTABK ( LEN , WRK2 , WRK1, WRK3 , 2) 

WRK4 (:,:,:) = BM (:,:,:, JM) + WRK3 (:,:,:) 
CALL INVERSE2 (LEN, WRK4, WRK1) 

CALL MULTAVK (LEN, WRK1, VRK1, VRK3 , 2) 
CALL MULTAVK (LEN, WRK1, VRK2 , VRK4 , 2) 

FN1 ( : , : , JM) = VRK3 ( : , :) 

FN2 ( : , : , JM) = VRK4 ( : , :) 

DO J = JM-1, 1, -1 

WRK1 ( : , : , : ) = ALPHA (:,:,:, J) 

VRK1 ( : , :) = FN1 ( : , :,J+1) 

VRK2 ( : , : ) = FN2 ( : , : , J+l) 

CALL MULTAVK (LEN, WRK1, VRK1, VRK3 , 2) 
CALL MULTAVK (LEN, WRK1 , VRK2 , VRK4 , 2) 

FN1 ( : , : , J) = VRK3 ( : , : ) + BETA1 < : , : , J) 
FN2 ( : , : , J ) = VRK4 ( : , :) + BETA2 ( : , :,J) 
ENDDO 


RETURN 

END 

SUBROUTINE SOMEAN ( DN1 , FN1, AM, BM, CM, JM) 
************************************************************* 

C Subroutine SOMEAN 

c ***..*ttt***t**tt**t*t*tt*#**«*****t*t*tt«*^«*ttt**^***t*t***t»* 

C WRITTEN BY: Yong Li , GSFC Code 910.3, Apr 1994 
C 

C PURPOSE: This routine solves the second-order ODE which originally results 
C the \del*4 equation on the sphere. But appplication is not so restricted 
C as far as the discretized equation is of the following format. 

C The unknown to be solved is a 2x1 vector (Lindzen and Kuo, 1969, MWR) . 

C This routine is used here to solve for the zonal mean, equivalent to FFT 
C coefficients (k = 0) in the \del A 4 direct solver. 

C 

C The discretized form of the ODE is: 

C 

C A ( j ) X(j-l) + B ( j ) X(j) + C(j) X(j + 1) = F ( j ) (1) 

C (for j = 2, 3, 4, JM-1) 

C 


C The boundary conditions are: 

C 

C A ( 1 ) X(l) + B ( 1 ) X ( 2 ) = F ( 1 ) (2.1) 

C 

C and A ( JM) X(JM-l) + B(JM) X(JM) = F(JM) (2.2) 

C 


C The coefficient A, B and C are all real. The forcing F and hence the 
C solution X also real. 

C 

C CALLED FROM: 

C DIRSOL4 

C 

C SUBROUTINES CALLED: 

C INVERSE: matrix inversion; 

C MULTAB: multiplication of two matrix; 

C MULTAV: multiplication of a matrix and a vector. 

C 

C VARIABLES: 


from 
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INPUT: 

AM-array storing the coefficient matrix A as in (1) ; 

BM-array storing the coefficient matrix B as in (1) ; 

CM-array storing the coefficient matrix C as in (1); 

DNl-array storing the vector of the forcing term F in (1) ; 
JM-number of total grid points including the two boundary points. 

OUTPUT: 

FNl-array storing the solution vector X of (1) ; 

LOCAL: 

ALPHA-array storing intermediate coefficient matrix alpha 
(Lindzen and Kuo, 1969); 

BETAl-real part of vector beta (Lindzen and Kuo, 1969); 

VRK1 -working array; 

VRK2 -working array; 

VRK3 -working array; 

VRK4 -working array; 

VRK5 -working array; 

VRK6 -working array; 

WRK1 -working array; 

WRK2 -working array; 

WRK3 -working array; 

WRK4 -working array; 

«****+**+*++*******+#+********+****+##****************#****+**********.* 
DIMENSION AM (2,2, JM) , BM(2,2,JM), CM(2,2,JM) 

#, DN1 ( 2 , JM) , FN1 ( 2 , JM) , ALPHA ( 2 , 2 , JM) , BETA1(2,JM) 

DIMENSION WRK1 (2,2) , WRK2(2,2), WRK3(2,2), WRK4(2,2) 

DIMENSION VRK1 (2) , VRK2(2), VRK3<2), VRK4<2), VRK5(2) 

#, VRK6 ( 2 ) 

WRK1 ( : , :) = AM ( : , : , 1 ) 

CALL INVERSEl (WRK1, WRK2) 

WRK1 ( : , : ) = - WRK2 ( : , : ) 

WRK4(:, :) = BM ( : , : , 1 ) 

CALL MULTAB (WRK1 , WRK4 , WRK3 , 2) 

ALPHA ( :,:,!) = WRK3 ( : , : ) 

VRK1 ( : ) = DN1 ( : , 1 ) 

CALL MULTAV ( WRK2 , VRK1, VRK3 , 2) 

BETA1 ( : , 1 ) = VRK3 ( : ) 

DO J = 2, JM- 1 

WRK1(:,:) = ALPHA ( : , : , J-l) 

WRK2( : , : ) = AM( : , : , J) 

CALL MULTAB (WRK2 , WRK1 , WRK3 , 2) 

WRK1 ( : , : ) = WRK3 ( : , : ) +BM(:,:,J) 

CALL INVERSEl (WRK1,WRK2) 

WRK1 ( : , :) = CM< : , :,J) 

CALL MULTAB ( WRK2 , WRK1 , WRK3 , 2) 

ALPHA ( : , : , J) = - WRK3 ( : , : ) 

WRK1 ( : , : ) = AM( : , : , J) 

VRK1 ( : ) = BETA1( : , J-l) 

CALL MULTAV (WRK1, VRK1, VRK3 , 2) 

VRK3 ( : ) = DN1 ( : , J) - VRK3 ( : ) 

CALL MULTAV (WRK2, VRK3 , VRK1, 2) 
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BETA1 ( : , J) = VRK1 ( : ) 
ENDDO 


VRK1 ( : ) = DN1 ( : , JM) 

VRK3 ( : ) = BETA1 { : , JM-1) 

WRK1 ( : , : ) = AM ( : , : , JM) 

CALL MULTAV ( WRK 1 , VRK3 f VRK5 , 2) 

VRK1 ( : ) = DN1 { : , JM) - VRK5 ( : ) 

WRK1(:,:) = ALPHA JM-1 ) 

WRK2 ( : , : ) = AM( : , : f JM) 

CALL MULTAB (WRK2 , WRK1 , WRK3 , 2) 

WRK4 ( : , : ) =BM(:,:,JM) + WRK3 < : , : ) 
CALL INVERSEl (WRK4 , WRK1) 

CALL MULTAV <WRK1, VRK1, VRK3 , 2) 


FN1(:,JM) = VRK3 ( : ) 

DO J = JM-1 , 1, -1 

WRK1 ( : , : ) = ALPHA ( : , : , J) 

VRK1 ( : ) = FN1 ( : , J+l) 

CALL MULTAV (WRK1 # VRK1, VRK3 , 2) 

FNl { : r J) = VRK3 ( : ) + BETA1(:,J) 
ENDDO 

RETURN 

END 

SUBROUTINE WRAP1 { IMA r IM r JM, A) 
DIMENSION A(IMA f JM) 

DO J = 1, JM 

DO I = IM+1, IMA 

A ( I , J) = A ( I - IM , J ) 

ENDDO 

ENDDO 

RETURN 

END 

SUBROUTINE WRAPD ( IMA , IM f JM , A) 
DIMENSION A ( IMA, JM) 

DO J = 1, JM 

A ( 1 , J) = A ( IM+ 1 , J ) 

DO I = IM+2, IMA 

A ( I # J ) = A ( I - IM , J ) 

ENDDO 

ENDDO 

RETURN 

END 
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Figure 2: (a) Spherical harmonic wave number one; (b) Same as (a), but with 2 Ax waves 
embedded; (c) After applying the semi-implicit formulation of fourth order diffusion to (b); 
(d) After applying the implicit formulation to (b); (e) After applying both the semi-implicit 
and the implicit formulations to (b). 
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Figure 2 (continued) 
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(e) 



Figure 2 (continued) 
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Figure 3: (a) Geopotential at 50 days for the integration with the initial condition of McDon- 
ald and Bates (1989) and a value of 0.03 for the uncentering parameter(contours every 60.0 
meters); (b) Same as (a) but with the fourth order diffusion applied to geopotential(contours 
every 60.0 meters); (c) Same as (a) but with the fourth order diffusion applied to divergence 
(contours every 60.0 meters); (d) Same as (a) but with the fourth order diffusion applied 
to both geopotential and divergence(contours every 60.0 meters). 
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Figure 3 (continued) 


(C) 



(d) 
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Figure 4: (a) Geopotential of the ECMWF analysis of 0000 UTC, 1 December, 1992 (con- 
tours every 60.0 meters); (b) Geopotential after 3 days of integration with a value of 0.05 
for the uncentering parameter(contours every 60.0 meters); (c) Same as (b) but with a value 
of 0.015 for the uncentering parameter and the the fourth order diffusion applied to both 
geopotential and divergence (contours every 60.0 meters). 
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(b) 



Figure 4 (continued) 
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Figure 5: 50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit 
GCM without the application of the horizontal diffusion. The value for the off-centering 
parameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c) 7-day forecast; (d) 10-day 
forecast. (Contour interval is 3 degrees). 





Figure 5 (continued) 












Figure 6: 50 mb temperature field from NASA/GLA 3-D semi-Lagrangian semi-implicit 
GCM with the application of the horizontal diffusion. The value for the off-centering pa- 
rameter is 0.05. (a) 3-day forecast; (b) 5-day forecast; (c) 7-day forecast; (d) 10-day forecast. 
(Contour interval is 3 degrees). 
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